Pair formation and collapse in imbalanced Fermion populations with unequal masses 
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We present an exact Quantum Monte Carlo study of the effect of unequal masses on pair formation 
in Fermionic systems with population imbalance loaded into optical lattices. We have considered 
three forms of the attractive interaction and find in all cases that the system is unstable and 
collapses as the mass difference increases and that the ground state becomes an inhomogeneous 
collapsed state. We also address the question of canonical vs grand canonical ensemble and its role, 
if any, in stabilizing certain phases. 
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I. INTRODUCTION 

The properties of fermionic systems with imbalanced 
populations have long been of interest for a variety of rea- 
sons. Of particular importance is the pairing mechanism 
resulting in composite bosonic particles and leading to su- 
perconductivity/superfluidity. The original focus 1 -^ was 
on polarized superconductors with an imbalance between 
the two spin populations where it was predicted that the 
Cooper pairs form with a center of mass momentum equal 
to the difference of the Fermi momenta of the two popu- 
lations. This Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
phased proved to be difficult to observe in solids and 
was only achieved recently in heavy fermion systems 3 -. 

The recent experimental realization of trapped ultra- 
cold fermionic atoms with tunable attractive interac- 
tions has intensified interest in this problem. Ex- 
periments, in which two hyperfine states of ultra-cold 
fermionic atoms play the role of up and down spins, 
have demonstrated the presence of pairing in the case 
of unequal populations* 4 ^ This has spurred much the- 
oretical work using a variety of approaches such as 
mean fie\S^^JMLlLl^±ll y effective Lagrangian^ 
bosonisationii and Bethe ansatz 1 ^ applied to the uniform 
system with extensions to the trapped system using the 
local density approximation (LDA). 

For the uniform case, the controversy continues on the 
details of the paired state: FFLO pairs forming with 
non-zero momentum vs breached pairing (BP) at zero 
momentu m 14 ' 15 ' 19 and on the robustness of such phases. 
Distinguishing features between these two possibilities in- 
clude the presence, for FFLO, of spatial modulations (in- 
homogeneity) in the pairing order parameter, and conse- 
quently a peak at non-zero momentum in the pair mo- 
mentum distribution function as opposed to the coexis- 
tence of a supcrfluid and normal component in a transla- 
tionally invariant and isotropic state in the BP scenario. 
In addition, excitations are gapped for the FFLO phase 
but not so for BP. 

No general consensus has emerged on which pairing 
type dominates. However, recent exact numerical work 
on the ground state of one-dimensional Fermi systems 



with imbalanced populations loaded into optical lattices 
has found that FFLO is the only pairing mechanism 
both in the unifor m 20 ' 21 ' 22 and the trappe d 20 ' 21 ' 23 ' 24 sys- 
tems. Similar results have been found without the opti- 
cal latticed On the other hand, approximate dynamical 
mean field (DMFT) results appear to indicate that no 
FFLO phase is present in the three-dimensional system^ 
In addition to the above equal-mass cases, of rele- 
vance in condensed matter and trapped ultra-cold atomic 
systems, there is great interest in the case of unequal 
masses. This arises naturally in trapped mixtures of dif- 
ferent atomic species, e.g. K-Rb, and also in cold dense 
quark matter where the c, b and t quarks are much heavier 
than the u, d and s quarks. It was argued theoretically 
that such mass and population imbalance leads to the 
BP phas o 27 ' 28 with gapless excitations which can be sta- 
bilized by longer-range interactions^ In this paper we 
examine this question with exact Quantum Monte Carlo 
(QMC) simulations. We focus on the one-dimensional 
system loaded into an optical lattice governed by the ex- 
tended Hubbard Hamiltonian 



W= - ^M4fi<T c z* + c L c j+i<r) 

la 

^ ^ Uij j( j(j' Tl{ a^lj a 1 ; 
i,j,a,a' 
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where Cj a (cj a ) are fermion creation (annihilation) oper- 
ators on spatial site j with the fermionic species labeled 



by g — 1,2 and 



l 3 " 



^ja^J a 



is the corresponding num- 



ber operator. The unequal masses are embodied in the 
unequal hopping parameters. We set the energy scale by 
taking t\ = 1 and studying the system as ti gets smaller 
(r»2 gets larger). In general, the interaction term Vij^a' 
can couple all fcrmions on all sites; in what follows we 
shall consider three different forms. 

For our simulations, we use a continuous imaginary 
time canonical "worm" (CW) algorithm where the total 
number of particles is maintained strictly constant^ In 
this algorithm, two worms are propagated, one for each 
type of fermion, which allows the calculation of the real- 
space Green functions of the two fermionic species, G a , 
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FIG. 1: (color online). The momentum distribution for the 
minority (left) and majority (right) populations. 

and the pair Green function, G pa i r , 

G pair (0 = <A. +i A}>, 

Aj = C j2 Cj!. (2) 

It is then immediate to obtain the momentum distribu- 
tions n a (k) and n pa i r (fc) from their Fourier transforms. 
We emphasize that this algorithm is exact: There are no 
approximations and the only errors are statistical (which 
typically correspond to the size of the symbols). The po- 
larization is given by P = (N 2 — Ni)/(N 2 + Ni), where 
Ni(N 2 ) is the minority (majority) particle numbers. The 
associated Fermi wave vectors are k-pa = TrN a /L. We ran 
our simulations at (3 = 2L, where L (typically L = 32) is 
the number of lattice sites. We have verified that these 
values arc large enough to yield the ground state proper- 
ties, which is our focus. Typical runs takes a day or two 
on a desktop computer. 

II. HEAVY MAJORITY: h > t 2 

We begin with the simplest form for the interaction 
term, Uij taa > = USij(l — 5 air7 >), i.e. only contact at- 
traction (U < 0) between the two species. Figure Q] 
shows what happens to the minority (left) and major- 
ity (right) momentum distributions as t 2 is decreased 
for fixed U = -3ii and N x = 9, N 2 = 11. For 
t 2 = 0.9ti, n a (k) drops sharply at the respective Fermi 
momenta, (fcp^, kp 2 ) while at the same time, the pair 
momentum distribution, n pa j r (fc), exhibits a maximum 
at k = \kp 1 — kp 2 \ (fig. [2] left) . This indicates pair for- 
mation at non-zero center-of-mass momentum and thus 
an FFLO phase i 20 ' 21 i 22 ' 23 i 24 As the mass of the major- 
ity particles increases, t 2 decreases, the attractive in- 
teraction effectively increases because \U\/t 2 increases. 
The increased binding is demonstrated by (nnni 2 ) which 
takes values between NiN 2 /L 2 (no binding) and N\/L 
(all minority particles have formed pairs )2£. This quan- 




k t 2 

FIG. 2: (color online). Left: The pair momentum distribution 
function. Right: The circles show the height of ripairdfcF! — 
kp 1) vs t2, the triangles show {nun^) scaled by 4 for better 
visibility. 



tity is shown in fig. [2] (triangles, right panel) scaled by 
4 for better visibility. Clearly, as t 2 decreases, (nnni 2 ) 
increases toward its upper limit. As a consequence of 
this increased binding, the FFLO effect first intensifies, 
its peak at k = \kp t — kp 2 \ increases in height, reaches a 
maximum for t 2 = 0.4^, then decreases and disappears 
as shown in fig. [5] (circles, right panel). Note, for exam- 
ple, that for t 2 = O.lti the peak of n pa ; r (fc), fig. [5] (left), 
is at k — 0. 

Also noteworthy is the appearance of a dip in n 2 (k), 
at k w kp 1 , whose depth increases with decreasing t 2 , 
reaching a maximum for t 2 = QAt\, corresponding to the 
maximum FFLO effect. Further decreasing t 2 washes 
this out. This feature can be understood as follows. As 
t 2 decreases, the interaction between the minority and 
majority effectively increases {\U\/t 2 ) 1 thus increasing the 
depletion of n a [k < kp a ). Most of the depletion for both 
ni(fc) and n 2 (k) happens for k < kp t since this is where 
most of the minority resides. This depletion for k < kp t 
leaves n 2 (k) with a bump for kp t < k < kp 2 . 

At first glance, it might seem that the disappearance of 
the FFLO peak, n pa ; r (fc = \kp t —kp 2 1), as t 2 decreases and 
the appearance of a peak n pa i r (/c = 0) signals pair forma- 
tion at k = and thus a BP phase. However, a simple 
argument offers another option. The Fermi momentum 
of the majority is kp 2 = nN 2 /L; consequently, as these 
particles get heavier, t 2 smaller, their kinetic energy be- 
comes negligible and can be ignored. To minimize its 
free energy, the system will then optimize the potential 
and kinetic energies of the light particles. The optimal 
potential energy is obtained when the light particles are 
on the same sites as the heavy ones. On the other hand, 
the kinetic energy is optimized when the light particles 
are delocalized. Both these energies can be optimized 
if the heavy particles coalesce, forming a contiguous re- 
gion with one heavy particle per site, jiq = 1. This 
region, then, acts as a platform on which the light par- 
ticles can delocalize over its entire extent while always 
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FIG. 3: (color online). Average density profiles for collapsed 
systems. Here U = —7ti and t-z = 0.04ii. Note the Friedel- 
like oscillations in nn. 
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FIG. 4: (color online). 8n versus t2/|i/| showing the collapse 
as t% gets small enough. 



being in contact with the heavy particles thus minimiz- 
ing their potential energy. This is indeed what happens 
as the density profiles, nn and n^, show in fig. [3] for two 
polarizations. To summarize, as t 2 decreases, the system 
goes from an FFLO phase to a spatially collapsed one. 

To quantify this collapse, we define the quantity 31 
(Sn) 2 = ( ( n u ) — Ni/L) 2 which is essentially zero when 
the system is uniform ((nn) = N±/L) and grows as the 
collapse happens. In fig. 2] we show 6n versus t2/\U\ 
for three polarizations and three couplings in each case. 
We see that, indeed, the system collapses as t 2 decreases 
and that for a given polarization, this collapse appears to 
happen at approximately the same value of t 2 /\U\. Also, 
the larger the polarization, the easier it is to trigger the 
collapse (larger t<x). This behavior holds for all the pa- 
rameters we examined, specifically —13 < U < — 1, and 
polarizations 0.1 < P < 0.55 

We, therefore, conclude that in the presence of only 
contact attraction, the BP phase is not realized as the 




FIG. 5: (color online) . Top: The pair momentum distribution 
for several values of t2 showing the disappearance of the FFLO 
peak. Bottom: Average density profile for same parameters 
as top panel but ti = 0.2ti showing the system when it first 
collapses. The contact interaction is U = — 4ti and the nn 
interaction is V = — 0.25t acting only between particles of 
opposite spin. 



majority population is made heavier; instead the system 
collapses. To counter the tendency of the heavy parti- 
cles to clump together as in fig. [3l we introduce longer- 
range interaction. It is reasonable to suppose that near 
neighbor (nn) repulsion between particles of the same 
species would tend to oppose such collapse. To this 
end we redid the above study but with the interaction 
term Hi = J2i Un n n l2 + V J2i a n i<y n i+ia with U < 
and V > 0. We studied this for 0.1 < P < 0.44 for 
— 10£i < U < —4ti and found that, indeed, the stability 
of the system against collapse is extended to smaller val- 
ues of *2 but that eventually the system always collapses. 
Furthermore, before the collapse, the system always ex- 
hibits FFLO pairing while after collapse, the nn repulsion 
term, V , produces density oscillations as the presence of 
near neighbors is opposed. Therefore, this second form 
of the interaction also fails to produce the BP phase. 

We now turn to the third interaction form which was 
proposed in Ref^ and argued to yield the BP phase. 
This form extends beyond the contact term and acts 
only between particles from different species. The idea 
is based on the assumption of three competing homoge- 
neous phases: a normal state of free fermions, a fully 
gapped BCS supcrfiuid and a gapless BP phase. It 
was argued that by giving the interaction term struc- 
ture in momentum space, one may be able to stabilize 
the BP phase; a Gaussian dependence on distance was 
used for the interaction potential^. In our simulations 
we cut the range at the nn level and, therefore, have 
Hi = J2i(Unnn.i2 + Vnnni + i2). Note that, unlike the 
previous form, V here acts only between n\ and n 2 . In 
momentum space, this has the form U + 2ycos(g) at mo- 
mentum q and, with U < and V < 0, is more attractive 
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FIG. 6: (color online). Top: 5n rises sharply at ii/|t/| « 0.03 
signalling a transition to an inhomogeneous density profile. 
Bottom: For ti/|(7| < 0.03, the system undergoes phase sep- 
aration as shown in this typical density profile. Part of the 
system is in a charge density wave phase while the remainder 
is in a free fermion phase. 




Minority, ^=-1.8 
Majority, n 2 =+0.72 
Pairs, grand can. 
Minority, N 1= 7 
Majority, N 2 =15 
Pairs, canonical 



L=32, U=-3t, P=20 



FIG. 7: (color online). Momentum distributions for the mi- 
nority, majority and pair populations using canonical and 
grand canonical QMC with ti — t-z. In the grand canonical 
case, (iVi) = 6.996 (/xi = -1.8), (N 2 ) = 14.949 (/x 2 = +0.72). 



III. HEAVY MINORITY: U < t 2 

So far, we have considered the case where the heavy 
particles are the majority which leads to collapsed config- 
urations like those in Fig. [3] The question then is: Will 
the system still collapse when the heavy particles are the 
minority and what form will the collapsed configurations 
take? We now consider this case with a heavy minority 
population Ni = 9, and a lighter majority, N2 = 13, on 
a 32-site lattice with j3 = 64. As before, we fix the light 
population hopping t% — 1 to define the energy scale and 
we study the collapse as a function of the heavy minority 
hopping parameter, ti/t^, and the attractive interaction, 
U jt-2 < 0. In the top panel of fig. [5] we show, like in fig.|H 
Sn as a function ti/|f/| for five different values of the in- 
teraction, U. The behavior is similar to that in fig. [U 
we find that for these values of Ni and N2, 5n increases 
sharply for <i/|C/| ps 0.03 signalling spatial collapse in 
the system. One candidate for the collapsed confiugra- 
tion is that, as before, the heavy particles (the minority 
in this case) form a contiguous region thus providing a 
platform for the lighter particles. This would then re- 
sult in a contiguous region with one heavy and one light 
particle per site and the excess light majority particles 
spread over the rest of the system. This, however, does 
not minimize the energy because the light particles resid- 
ing on the heavy particle platform are blocked: They are 
in a Mott state and have zero kinetic energy. In the bot- 
tom panel of fig. [5] we show a typical density profile of a 
collapsed configuration. It is easy to understand energet- 
ically why this density wave structure is favored over the 
previous candidate: The light particles are never blocked 
in a Mott region and can always hop to neighboring sites 
to optimize the free energy. 

We note that the configuration in fig. [5] corresponds 
to phase separation: In one region the system is in a 
charge density wave phase, in the other region it has free 
fermions. These two phases co-exist. On the other hand, 
the configurations in fig.[3]correpsond to spatial collapse: 
The system has a regions void of particles and regions 
where all the particles reside. 



IV. CANONICAL VS GRAND CANONICAL 



at small q in contradistinction to the pure contact inter- 
action which is independent of q. Our exact QMC results 
show, however, that the system is still unstable and col- 
lapses as ti decreases, just like in the contact interaction 
case. This is shown in fig. [5j No BP phase has been 
found. In fact, with the longer-range attraction, the sys- 
tem is even more unstable as can be easily understood 
by the same real space argument presented above for the 
contact interaction case. 



Finally, we consider the question of canonical versus 
grand canonical ensembles. It has been suggeste d 15 i 29 
that the stability of BCS, FFLO or BP phases can de- 
pend on whether one fixes the populations or the chem- 
ical potentials. All the above results have been ob- 
tained with the CW algorithm where N\ and N2 are 
kept strictly fixed. An alternative is to add the chem- 
ical potential term, Mi^ii — 1^2^2)1 to the Hamil- 
tonian, eq.fl]), and use a grand canonical QMC algo- 
rithm such as the Determinant Quantum Monte Carlo 
(DQMC) algorithm^ We compare in fig. [7] the momen- 
tum distributions obtained with the CW and the DQMC 
algorithms with \x\ and /Z2 chosen to give average fill- 



5 



ings corresponding to the fixed fillings in the canonical 
case. There is no disagreement between the two; in par- 
ticular, the FFLO phase is seen to be stable whether 
one fixes the populations or the chemical potentials. Of 
course, if one changes ti while holding the chemical po- 
tentials fixed, the average populations, (Ni) and (N2), 
will change. Nonetheless, whatever (N\) and (N2) one 
has obtained by fixing fi\ and 112, one will obtain the same 
physics in the canonical ensemble by fixing the popula- 
tions to the corresponding values, as seen in fig. [7] One 
is free to study these phases and their stability in either 
ensemble. 



V. CONCLUSIONS 

In summary, we have studied, using exact QMC simu- 
lations, the effect of mass differences between two unbal- 
anced Fermion populations. For the case where the ma- 
jority population is the heavier, we performed our study 
with three possible attractive interaction terms. In all 
three cases, we have found the system to be unstable and 



to collapse when the mass disparity is large enough: The 
BP phase is not realized by tuning the mass ratio between 
the two populations. For the case where the minority is 
heavier, we showed the system still collapses when the 
mass difference is large enough, but in this case it forms 
density wave structures. We have also shown that fixing 
the populations or the chemical potentials leads to the 
same physics. The stabilization of sought-after phases is 
not favored by one ensemble rather than another. 
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